I Miloncini¶
| Last Name | First Name | Student Number |
|---|---|---|
| Bertinelli | Gabriele | 2103359 |
| Bonato | Diego | 2091250 |
| Di Prima | Giacomo | 2086992 |
| Viterbo | Giuseppe | 2086516 |
Part 2 - Data Analysis¶
The N-body Problem¶
The N-body problem refers to the challenge of predicting the movement of a collection of celestial bodies that are gravitationally interacting. This issue holds significant relevance in the field of astrophysics due to its wide application across various astrophysical contexts and scales.
The N-body problem is an extension of Newton's law of universal gravitation problem.
The acceleration of the body $i$ due to all the other bodies $j$ is expressed by $$ \begin{equation*} \frac{ d^2 \vec{x}_i } {dt^2} = - G \sum^N_{j=1, j\neq i} m_j \frac{ \vec{x}_i - \vec{x}_j }{ | \vec{x}_i - \vec{x}_j |^3 } . \end{equation*} $$
Each N-body system possesses six constants of motion, comprising the position and velocity of the center of mass, along with four integrals of motion encompassing energy and the three components of angular momentum.
While an analytic solution exists for N=2, extending it to a general N-body system remains elusive. Adding just a third particle increases the unknowns to 18, necessitating constraints for an analytic solution, which is only feasible for specific cases like the circular restricted three-body problem.
Numerical methods are indispensable for integrating the system's differential equations. Notably, the problem's numerical complexity scales at O(N^2), posing a significant hurdle in system modeling due to rapid growth with increasing N.
N-body Units¶
N-body units offer a practical framework for N-body simulations, simplifying calculations by setting $G=1$. Conversion to physical units can be done retroactively, incorporating typical values relevant to the astrophysical system under study. However, this approach becomes inadequate when incorporating data on stars' physical parameters, stellar evolution, supernova explosions, etc. Nonetheless, treating stars as point masses allows for exploiting the scale invariance of N-body simulations.
Henon Units¶
Additionally, Henon Units provide another set of normalized units used in N-body simulations.
In Henon Units, several fundamental physical quantities are set to unity, simplifying the equations of motion and standardizing the numerical treatment of N-body systems.
$G=1$: this simplifies the gravitational force equation, as the force between two masses is directly proportional to their product.
The total mass of the N-body system is also set to 1: this normalization simplifies mass-related calculations and eliminates the need to account for varying masses across different bodies.
With this scaling the total energy of the system is set to $-0.25$: this normalization allows for straightforward comparison of energy-related quantities and ensures consistency across simulations. Specifically, the kinetic energy is set to $0.25$ and the potential energy to $-0.5$.
By normalizing these fundamental quantities the equations of motion for the N-body system become simpler, facilitating numerical integration and analysis. With all simulations using the same normalized units, it becomes easier to compare results across different studies and simulations.
In our analysis we worked in these units, anyway, physical units can be easily recovered with the functions included in the module fireworks.nunits.
Project Description¶
Stellar clusters orbiting in and around a galaxy are subject to the tidal field of that galaxy. Comparing the potential field generated by the galaxy and the one generated by the stellar cluster it is possible to define a tidal radius $$ \begin{equation} r_t \approx r_G \biggl(\frac{m_{cl}}{M_G}\biggr)^{\frac{1}{3}} \end{equation} $$ where $r_G$ is the distance of the stellar cluster from the center of the galaxy, $M_G$ is the mass of the galaxy within $r_G$ and $m_{cl}$ is the mass of the stellar cluster. At $r > 2r_t$ the potential generated by the galaxy starts to dominate over the potential of the cluster, therefore we can consider the star escaped by the stellar system. Therefore, along the orbit, the stellar cluster loses mass creating stellar streams. In some cases, the stellar cluster is destroyed.
Plummer Sphere¶
In this work, we run different simulations of a stellar cluster losing mass while orbiting around a point-mass external potential.
The initial conditions have been drawn from a self-gravitating Plummer Profile in virial equilibrium composed of a different number of stars. All units are Henon Units.
We tested different configurations of the system:
- 500, 1000, 5000 stars
- eccentricity: 0.0, 0.5, 0.9
We then decided to analyze and present in this report the realization with 1000 stars and eccentricity equal to 0.0 and 0.5.
Point-mass Potential¶
We included in the simulation a simple external potential available in Fireworks: fireworks.nbodylib.potentials.Point_Mass. We set the mass of the point-mass to $10^3$, in this way, the total mass of the external potential is $10^3$ times larger than the stellar cluster mass.
The Fireworks Python Package¶
Fireworks is a Python library that contains tools to initialize and evolve N-body systems and it can be used to simulate collisionless systems, collisional systems and orbit integration. It also contains some other useful tools as pyfalcon and TSUNAMI.
Fireworks contains the following submodules:
particles: it contains the classParticlesin which all the information about the particles can be stored (position, velocity, mass), and it provides useful tools to retrieve physical quantities of the system (e.g. CoM quantities, potential, kinetic and total energies);ic: it contains functions through which initial conditions can be generated. The module contains two functions:ic_random_normal, which draws initial conditions from a normal distribution, andic_two_body, which generates initial conditions for the case of a two-body system;Nbodylib: it contains functions to estimate the gravitational forces and accelerations. It is divided into other four modules:dynamics: it contains functions to estimate accelerations due to gravitational forces.acceleration_pyfalcon, used in our work, is a function of this module. In this module we can find our ad hoc built acceleration estimate functions:acceleration_directandacceleration_direct_vectorized. The first computes gravitational acceleration between particles using a direct method (i.e. for loops), while the latter leverages the broadcasting operations ofnumpy.array.integrators: it contains integrators used to integrate the ODE equations of the motion and evolve the system in time.integrator_leapfrog, used in our work, is a function of this module;nunits: it contains the classNbody_units, necessary to transform the data from physical units to N-body units or vice-versa;potentials: it contains a collection of classes and functions to estimate acceleration due to gravitational forces of a fixed potential. In this module, we can find some classes among themPotential_Base, used to initialize new potentials, andPoint_Mass, which assumes the presence of a point of mass M fixed at the center of the frame of reference.timesteps: it contains functions to estimate the adaptive timestep for the N-body integrations.
Simulations¶
The ad-hoc function built for the project can be found in the file simulation.ipynb. Here we present some key points.
First, we had to put the stellar clusters in orbit around the external potential. After having moved the stars to the frame of reference of the center of mass (CoM) of the cluster, we moved the stars to the initial position
particles.pos[:, 0] += initial_position * np.ones(len(particles))
so that the center of mass coordinates were $x=$ initial_position/2 and $y=z=0$.
Then we used the function ic_two_body to calculate the velocity of the cluster based on the eccentricity of the orbit, updating the velocity of the stars along the y-axis
particles.vel[:, 1] += vel_cluster
We implemented the function softening to calculate the softening length as the mean distance between the stars. We passed this value to the function that generates the point-mass potential
potential_PointMass = fnp.Point_Mass(Mass=M_G, softening = mean_distance)
We choose as parameters for the simulation:
initial_position = 10 # Initial position along x direction
tstep = 0.01 # Time step for integration
N_orbit = 2 # (or 3) Number of orbits to integrate
Analysis¶
The code used to produce the plots below can be found in the file Analysis.ipynb.
Here below we present a description and an analysis of the results found for two different systems:
- Number of stars = 1000, initial position = 10 (along x direction)
- eccentricity = 0.0
- eccentricity = 0.5
a) N=1000, e=0.0¶
Figure a.1. This plot shows the distance between the "true" center of mass (CoM) $r_{CoM}$ of the particles and the iterative one $r_{CoM, tidal}$. The first is the CoM found taking the mean of the position of each particle. The latter is found iteratively in the following way: use the stars considered bounded in the previous snapshot to calculate the CoM, update the tidal radius in the new position of the center of mass and consider as bounded stars only those within 2 tidal radii, $2\, r_t$, and repeat this process till convergence.
We see that the two centers of mass deviate greatly as the simulation time passes. This is because when the mass loss starts to be significant the center of mass can be not centered in the density peak of the system.
Figure a.2. In this plot, a projection of the orbit on the X-Y plane (left) and on the Y-Z plane (right) is presented. In the center (i.e. in $(0,0)$), indicated by the green cross, there is the origin of Point Mass (PM) potential.
On the left side, we can see how the "true" CoM (small star) and the iterative CoM (small circle) deviate along the orbit and as the simulation time passes, indicated by the colorbar. The integrated orbit, relative to the iterative CoM, is pretty stable on the X-Y plane.
On the right side, some small fluctuations are visible, but the orbit can anyway be considered as planar.
From now on, the iterative CoM will be referred as CoM.
Figure a.3. On the left, we can see the mass loss at each time step. On the right, is the cumulative mass loss.
We can see that the cluster loses almost all the mass ($>80\%$) within 500 time steps, due to tidal effects.
Figure a.4. In this plot, the radii within which 10%, 50% and 90% of the mass is contained are displayed. For this plot we considered the mass contained inside $2\, r_t$, at each snapshot.
It is clear that until time step 380, the internal region $r<$ 10% mass radius and the more external ones are decoupled. This means that while the external regions expand due to tidal forces and lose mass (making the 90% mass radius shrink as time passes), the internal region remains roughly the same size.
Figure a.5. In this plot, we show the density profile inside 5 radial shells of volume. For the mass, we considered only the stars within $2\, r_t$. We show different time steps using different colors.
As time passes, the peak of the density profile stays pretty much at the same radii position, until the major mass loss event at $>300$ time steps occurs. The peak shifts towards smaller radius values because a lot of mass is lost from the external regions and thus the density peaks end to coincide with the 10% mass distribution, ~ 0.18 $|r-r_{CoM}|_{shell}$.
Figure a.6. In this figure, we plot the tangential velocity dispersion w.r.t. the CoM, contained in the corresponding half mass radius. We show In different colors the integration time steps.
We can see in the top right corner, before 300 time steps, three distinct phases:
- In the first one, until ~ 100 time steps, the dispersion grows very rapidly in the inner region;
- In the second phase, until ~250 time steps, the inner region starts to become more homogeneous in the tangential direction, because it expands due to the tidal forces, but the stars are still within $2\, r_t$.;
- In the third one, until ~400 time steps, the half mass radius shrinks due to mass loss while velocity dispersion remains almost the same because only the outer layers are stripped.
Figure a.7. Here we show some snapshots of the evolution of the system on the X-Y plane. We highlight in orange the stars within $2\, r_t$ at each snapshot. We see that when the massive mass loss occurs, at ~400 time steps, and further the stars distribute along the orbit becoming the actual stellar stream. Towards the end of the evolution, the stars fall in a smaller on average orbit around the PM.
Figure a.8. Here we show some snapshots of the evolution of the system on the X-Y plane, selecting only the stars belonging to the tidal stream (i.e. $r>2\,r_t$). We highlight with colors the number density of the stars, found with a 2D kernel density estimation. In red the region with the highest density. Until before the total disruption of the cluster, before ~500 time steps, the density peaks are closer to the CoM because that is where the stripped stars are coming from. In the final snapshot, when the cluster is already disrupted, the density is almost homogeneous along the tidal streams.
Figure a.9. Here we show some snapshots of the evolution of the system on the X-Y plane. We highlight with colors the differential tangential velocity of each star, w.r.t the CoM. Due to the tidal forces, stars from the inner and outer parts, from the center of the cluster, leave the cluster with faster and slower velocities w.r.t. the CoM respectively. Thus, the inner stars become a leading arm of the stellar stream while the outer stars form a trailing arm.
Higher velocities are shown in bluer colors, while lower ones are in redder colors, highlighting the structure of the stellar stream.
Figure a.10. This plot is similar to the one above (Fig. 9). There are two main differences:
- In the left plots, we show only stars in the stellar stream (i.e. those having $r>2\, r_t$).
- On the right, we describe the distribution of the velocities of the two arms (Leading in blue, Trailing in red), comparing them with the velocities distribution of the stars inside $2\, r_t$.
In order to separate the arms, we used the following selection rule:
- Leading stream: $\phi > \phi_{CoM}$
and$\rho < \rho_{CoM}$ - Trailing stream: $\phi < \phi_{CoM}$
and$\rho > \rho_{CoM}$
The reference system for the cylindrical coordinates $(\rho, \phi)$ is centered in $(X,Y)=(0,0)$. $\phi$ grows anti-clockwise, while $\rho$ is the usual radial distance from the center of the reference frame.
Figure a.11.
Top panel: Here we show some snapshots of the evolution of the system on the X-Y plane. We plotted all the stars in the system.
The solid (ticked) black line shows the region inside $2\,r_t$ ($r_t$). The colorbar indicates the total mutual energy between stars inside $2\,r_t$ (without considering the external PM potential). The total energy is defined as $E_{tot}=E_{kin}+E_{pot, mutual}$.
In the legend, we added two estimates of the percentual errors, defined as follows
$$ \begin{equation*} Err_{\%} = \frac{N_{R<(\alpha\cdot\,r_t)} - N_{E_{tot}<0}}{N_{E_{tot}<0}} \quad , \quad \alpha \in \{1, 2\} \end{equation*} $$With $N_{E_{tot}<0}$ the number of all the stars with a negative total energy and $N_{R<(\alpha\cdot\,r_t)}$ the number of stars inside 1 (or 2) $r_t$.
A negative value of $Err_{\%}$ means we are underestimating the number of stars still dynamically bound to the cluster.
Botto panel:
As we can see from the histogram, one $r_t$, on average, underestimates the number of stars still bound to the cluster. While $2\,r_t$ both underestimates and overestimates the number of stars. So we can conclude that the conservative condition of $2\,r_t$, in a simulation setup where energies are available, is not proper in defining the limit inside which a cluster is gravitationally bounded. A more precise limit is given by the particles' total energy.
b) N=1000, e=0.5¶
In this Section, we propose a comparison of the results with the ones discussed above. Therefore the plots will be not described again.
Figure b.1. As is Fig. a.1 we can see how the "true" CoM deviates from the iterative CoM.
Figure b.2. The system stays stable and planar for all the evolution time. The iterative CoM traces the correct evolution of the cluster along the orbit, around the PM.
Figure b.3. In this system the cumulative mass loss is less steep because the cluster loses mass in a longer time range due to the fact that is along a more eccentric orbit.
Figure b.4. We can see, similarly to the case above, that the inner layers are decoupled from the outer layers until the cumulative mass loss starts to be significant (i.e. above 1400 time steps).
We are looking at the same phenomenon, as in the case $e=0.0$, but on a longer timescale.
Figure b.5. In this case, the density peak moves towards outer shells. This effect is because the stars in the cluster move toward outer shells but are not stripped from the host body, as fast as in the case of $e=0.0$ (i.e. the cluster retains its mass for longer time).
Figure b.6. We see the same situation as the case $e=0.0$ until time step 400. It's the same result just on another timescale.
Figure b.7. Here we show some snapshots of the evolution of the system on the X-Y plane. We highlight in orange the stars within $2\, r_t$ at each snapshot. Thanks to the eccentric orbit, the cluster retains its mass for a longer time, surviving almost until the end of the simulation (after time step ~1400, almost all the mass is lost).
Figure b.8. Here we show some snapshots of the evolution of the system on the X-Y plane, selecting only the stars belonging to the tidal stream (i.e. $r>2\,r_t$). We highlight with colors the number density of the stars, found with a 2D kernel density estimation.
As said above, the cluster retains more mass almost until the end of the simulation.
Figure b.9. There are no meaningful differences in the results w.r.t. the case above, except for the fact that the cluster remains visible almost until the end of the simulation.
Figure b.10. From the plots above, we can appreciate the fact that leading and trailing differential tangential velocity distributions start to separate significantly from the cluster distribution only towards the end of the simulation.
A statistical analysis should be run in order to study the overlapping of the distributions, considering, for example, the 95% credibility interval.
Figure b.11. In the case of an eccentric orbit we see that the limit $r_t$ underestimates the number of stars, while $2\,r_t$ overestimates it. So an intermediated radius
$r_t<r<2\,r_t$ can be found in order to estimate which stars are bound to a cluster.
Conclusions¶
From the comparison of the two cases we can draw several conclusions:
Eccentric systems can retain more mass for longer times;
For circular orbits, the conservative condition of $2\,r_t$ gives a lower bound for the gravitationally bound condition. A more precise limit is given by the particles' total energy.
While for eccentric orbits, $r_t$ and $2\, r_t$ give a lower and an upper bound respectively, so an intermediated radius $r_t<r<2\,r_t$ can be found as a proxy for the gravitationally bound condition.In both cases, the inner layers are decoupled from the external ones. This means that while the external regions expand due to tidal forces and lose mass (making the 90% mass radius shrink as time passes), the internal regions remain roughly the same size. The process is the same but with different timescales, depending on the moment the cumulative mass loss becomes to be significant (i.e. $>80\%$).
The dispersion of the tangential velocity, in Fig. 6, presents a similar behavior in both cases, only on different timescales. The three phases, that we identified, are a fast increase in velocity dispersion, a mixing of the inner layers' velocities and a stripping of the outer layers.
Future works should investigate the dependencies of using the tidal radius condition for electing gravitationally bound stars, depending on the eccentricities of the orbits.
A more robust statistical analysis should be run in order to study the overlapping of the tangential velocities distribution, considering, for example, the 95% credibility interval.